In [1]:
%matplotlib inline

In [2]:
from ndreg import *
import matplotlib
import ndio.remote.neurodata as neurodata

In [13]:
import clarityviz as cl
import numpy as np
import nibabel as nib

In [4]:
reload(cl)


Out[4]:
<module 'clarityviz' from 'clarityviz/__init__.pyc'>

In [15]:
path = 's275_to_ara3_regis.nii'

In [21]:
im = nib.load(path)
im = im.get_data()
print(im.shape)


(1140, 800, 1320)

In [36]:
im_slice = im[:,:,660]

plt.imshow(im_slice, cmap='gray')
plt.show()



In [32]:
im_slice = im[600,201:300,201:400]

plt.imshow(im_slice, cmap='gray')
plt.show()



In [87]:
im_slice = im[600,201:300,201:400]

plt.imshow(im_slice, cmap='gray')
plt.show()


Extracting the top 10% brightest points


In [105]:
points = extract(im_slice, b_percentile = 0.9)


20245.29
(This will take couple minutes)
Above threshold=1964
total points: 1964
Output num points: 1964
Finished

In [111]:
print(im_slice.shape)


(99, 199)

In [106]:
print(points)


[[   16   198 22334]
 [   31    99 20485]
 [   35    96 22033]
 ..., 
 [   98   196 21508]
 [   98   197 25175]
 [   98   198 26081]]

In [109]:
hist,bins = np.histogram(im_slice.flatten())

In [110]:
print(max(bins))


28510.0

In [113]:
bim = np.zeros(im_slice.shape)

In [114]:
for point in points:
    bim[point[0], point[1]] = point[2]

In [115]:
plt.imshow(bim, cmap='gray')
plt.show()



In [116]:
plt.imshow(im_slice, cmap='gray')
plt.show()



In [71]:


In [41]:
def extract_bright_points(im, num_points=10000, b_percentile=0.75, optimize=True):
    """
    Function to extract points data from a np array representing a brain (i.e. an object loaded
    from a .nii file).
    :param im: The image array.
    :param num_points: The desired number of points to be downsampled to.
    :param b_percentile: The brightness percentile.
    :param optimize:
    :return: The bright points in a np array.
    """
    # obtaining threshold
    (values, bins) = np.histogram(im, bins=1000)
    cumValues = np.cumsum(values).astype(float)
    cumValues = (cumValues - cumValues.min()) / cumValues.ptp()

    maxIndex = np.argmax(cumValues > b_percentile) - 1
    threshold = bins[maxIndex]
    print(threshold)

    total = im.shape[0] * im.shape[1] * im.shape[2]
    #     print("Coverting to points...\ntoken=%s\ntotal=%d\nmax=%f\nthreshold=%f\nnum_points=%d" \
    #           %(self._token,total,self._max,threshold,num_points))
    print("(This will take couple minutes)")
    # threshold
    im_max = np.max(im)
    filt = im > threshold
    # a is just a container to hold another value for ValueError: too many values to unpack
    # x, y, z, a = np.where(filt)
    t = np.where(filt)
    x = t[0]
    y = t[1]
    z = t[2]
    v = im[filt]
    #     if optimize:
    #         self.discardImg()
    #     v = np.int16(255 * (np.float32(v) / np.float32(self._max)))
    l = v.shape
    print("Above threshold=%d" % (l))
    # sample

    total_points = l[0]
    print('total points: %d' % total_points)

    if not 0 <= num_points <= total_points:
        raise ValueError("Number of points given should be at most equal to total points: %d" % total_points)
    fraction = num_points / float(total_points)

    if fraction < 1.0:
        # np.random.random returns random floats in the half-open interval [0.0, 1.0)
        filt = np.random.random(size=l) < fraction
        print('v.shape:')
        print(l)
        #         print('x.size before downsample: %d' % x.size)
        #         print('y.size before downsample: %d' % y.size)
        #         print('z.size before downsample: %d' % z.size)
        print('v.size before downsample: %d' % v.size)
        x = x[filt]
        y = y[filt]
        z = z[filt]
        v = v[filt]
        #         print('x.size after downsample: %d' % x.size)
        #         print('y.size after downsample: %d' % y.size)
        #         print('z.size after downsample: %d' % z.size)
        print('v.size after downsample: %d' % v.size)
    points = np.vstack([x, y, z, v])
    points = np.transpose(points)
    print("Output num points: %d" % (points.shape[0]))
    print("Finished")
    return points

In [90]:
def extract(im, num_points=10000, b_percentile=0.75, optimize=True):
    """
    Function to extract points data from a np array representing a brain (i.e. an object loaded
    from a .nii file).
    :param im: The image array.
    :param num_points: The desired number of points to be downsampled to.
    :param b_percentile: The brightness percentile.
    :param optimize:
    :return: The bright points in a np array.
    """
    # obtaining threshold
    (values, bins) = np.histogram(im, bins=1000)
    cumValues = np.cumsum(values).astype(float)
    cumValues = (cumValues - cumValues.min()) / cumValues.ptp()

    maxIndex = np.argmax(cumValues > b_percentile) - 1
    threshold = bins[maxIndex]
    print(threshold)

#     total = im.shape[0] * im.shape[1] * im.shape[2]
    total = im.shape[0] * im.shape[1]

    #     print("Coverting to points...\ntoken=%s\ntotal=%d\nmax=%f\nthreshold=%f\nnum_points=%d" \
    #           %(self._token,total,self._max,threshold,num_points))
    print("(This will take couple minutes)")
    # threshold
    im_max = np.max(im)
    filt = im > threshold
    # a is just a container to hold another value for ValueError: too many values to unpack
    # x, y, z, a = np.where(filt)
    t = np.where(filt)
    x = t[0]
    y = t[1]
#     z = t[2]
    v = im[filt]
    #     if optimize:
    #         self.discardImg()
    #     v = np.int16(255 * (np.float32(v) / np.float32(self._max)))
    l = v.shape
    print("Above threshold=%d" % (l))
    # sample

    total_points = l[0]
    print('total points: %d' % total_points)

#     if not 0 <= num_points <= total_points:
#         raise ValueError("Number of points given should be at most equal to total points: %d" % total_points)
#     fraction = num_points / float(total_points)

#     if fraction < 1.0:
#         # np.random.random returns random floats in the half-open interval [0.0, 1.0)
#         filt = np.random.random(size=l) < fraction
#         print('v.shape:')
#         print(l)
#         #         print('x.size before downsample: %d' % x.size)
#         #         print('y.size before downsample: %d' % y.size)
#         #         print('z.size before downsample: %d' % z.size)
#         print('v.size before downsample: %d' % v.size)
#         x = x[filt]
#         y = y[filt]
#         z = z[filt]
#         v = v[filt]
#         #         print('x.size after downsample: %d' % x.size)
#         #         print('y.size after downsample: %d' % y.size)
#         #         print('z.size after downsample: %d' % z.size)
#         print('v.size after downsample: %d' % v.size)
    points = np.vstack([x, y, v])
#     points = np.vstack([x, y, z, v])
    points = np.transpose(points)
    print("Output num points: %d" % (points.shape[0]))
    print("Finished")
    return points

In [17]:
def plot_hist(im, title=''):
    hist,bins = np.histogram(im.flatten(),256,[0,256])

    cdf = hist.cumsum()
    cdf_normalized = cdf * hist.max()/ cdf.max()

    plt.plot(cdf_normalized, color = 'b')
    plt.hist(im.flatten(),256,[0,256], color = 'r')
    plt.title(title)
    plt.xlim([0,256])
    plt.legend(('cdf','histogram'), loc = 'upper left')
    plt.show()

In [ ]:


In [33]:
fig = plt.figure()
a=fig.add_subplot(1,2,1)
imgplot = plt.imshow(im[:,:,1000])
a.set_title('Before')
plt.colorbar(ticks=[0.1,0.3,0.5,0.7], orientation ='horizontal')
a=fig.add_subplot(1,2,2)
imgplot = plt.imshow(im[600,:,:])
imgplot.set_clim(0.0,0.7)
a.set_title('After')
plt.colorbar(ticks=[0.1,0.3,0.5,0.7], orientation='horizontal')


Out[33]:
<matplotlib.colorbar.Colorbar at 0x7f911e922d50>